#include "fermiqcd.h"

void main(int argc, char** argv) {
  string info=argv[1];
  mpi.open_wormholes(argc, argv);
  define_base_matrices("FERMILAB");

  FILE *fp;
  int box[4];
  coefficients coeff_light, coeff_heavy;

  string gauge_filename=          prompt(info, "INPUT_GAUGE_FILENAME", "COLD");
  string fermi_filename=          prompt(info, "OUTPUT_FERMI_PREFIX", "fermi");

  box[0]=               (int) val(prompt(info, "T", "16"));
  box[1]=               (int) val(prompt(info, "X", "4"));
  box[2]=               (int) val(prompt(info, "Y", "4"));
  box[3]=               (int) val(prompt(info, "Z", "4"));

  string action_light=          prompt(info, "LACTION", "CLOVER");
  string inverter_light=        prompt(info, "LINVERTER", "MINRES");
  mdp_real aprecision_light=val(prompt(info, "LABSOLUTE_PRECISION", "0"));
  mdp_real rprecision_light=val(prompt(info, "LRELATIVE_PRECISION", "1e-6"));

  coeff_light["kappa_s"]= val(prompt(info, "LKAPPA_S", "0.13"));
  coeff_light["kappa_t"]= val(prompt(info, "LKAPPA_T", "0.13"));
  coeff_light["r_s"]=     val(prompt(info, "LR_S", "1.0"));
  coeff_light["r_t"]=     val(prompt(info, "LR_T", "1.0")); // this must me 1  for sse
  coeff_light["c_{sw}"]=  val(prompt(info, "LC_SW", "0.0"));
  coeff_light["c_E"]=     val(prompt(info, "LC_E", "1.0"));
  coeff_light["c_B"]=     val(prompt(info, "LC_B", "1.0"));

  string action_heavy=          prompt(info, "HACTION", "CLOVER");
  string inverter_heavy=        prompt(info, "HINVERTER", "MINRES");
  mdp_real aprecision_heavy=val(prompt(info, "HABSOLUTE_PRECISION", "0"));
  mdp_real rprecision_heavy=val(prompt(info, "HRELATIVE_PRECISION", "1e-6"));

  coeff_heavy["kappa_s"]= val(prompt(info, "HKAPPA_S", "0.13"));
  coeff_heavy["kappa_t"]= val(prompt(info, "HKAPPA_T", "0.13"));
  coeff_heavy["r_s"]=     val(prompt(info, "HR_S", "1.0"));
  coeff_heavy["r_t"]=     val(prompt(info, "HR_T", "1.0")); // this must me 1  for sse
  coeff_heavy["c_{sw}"]=  val(prompt(info, "HC_SW", "0.0"));
  coeff_heavy["c_E"]=     val(prompt(info, "HC_E", "1.0"));
  coeff_heavy["c_B"]=     val(prompt(info, "HC_B", "1.0"));

  string smearing_type=prompt(info, "SMEARING", "NONE");
		 
  char tmp[100];
  int i,j,k,a,b,c,t;
  mdp_lattice lattice(4,box,default_partitioning0,torus_topology,0,1,false);
  gauge_field U(lattice,3);
  fermi_field psi(lattice,3);
  fermi_field chi(lattice,3);
  fermi_field phi(lattice,3);
  site x(lattice);
  mdp_array<double,1> sum(box[0]);

  if(gauge_filename=="COLD") set_cold(U);
  else U.load(gauge_filename);
  if(coeff_light.has_key("c_{sw}") ||
     coeff_heavy.has_key("c_{sw}")) compute_em_field(U);
  
  for(a=0; a<psi.nspin; a++)
    for(i=0; i<psi.nc; i++) 

      for(b=0; b<psi.nspin; b++)
	for(j=0; j<psi.nc; j++) {

	  if(action_light=="CLOVER") 
	    default_fermi_action=FermiCloverActionFast::mul_Q;
#if defined(SSE2)
	  else if(action_light=="CLOVERSSE") 
	    default_fermi_action=FermiCloverActionSSE2::mul_Q;
#endif
	  else
	    error("unknown action");
	  
	  if(inverter_light=="MINRES") 
	    default_fermi_inverter=MinimumResidueInverter<fermi_field,gauge_field>;
	  else if(inverter_light=="BICGSTAB") 
	    default_fermi_inverter=BiConjugateGradientStabilizedInverter<fermi_field,gauge_field>;
	  else
	    error("unknown inverter");


	  sprintf(tmp,"%s_%f_%s_%i%i.mdp", gauge_filename.c_str(),
		  coeff_light["kappa_s"], smearing_type.c_str(),a,i);
	  
	  if(is_file(tmp)) chi.load(tmp);
	  else {
	    psi=0;
	    if(on_which_process(lattice,0,0,0,0)==ME) { x.set(0); psi(x,a,i)=1; }
	    psi.update();
	    inversion_stats stats=
	      mul_invQ(chi,psi,U,coeff_light,
		       aprecision_light,
		       rprecision_light);
	    if(smearing_type=="NONE") mpi << "Local smearing\n";

	    chi.save(tmp);
	  }


	  if(action_heavy=="CLOVER") 
	    default_fermi_action=FermiCloverActionFast::mul_Q;
#if defined(SSE2)
	  else if(action_heavy=="CLOVERSSE") 
	    default_fermi_action=FermiCloverActionSSE2::mul_Q;
#endif
	  else
	    error("unknown action");
	  
	  if(inverter_heavy=="MINRES") 
	    default_fermi_inverter=MinimumResidueInverter<fermi_field,gauge_field>;
	  else if(inverter_heavy=="BICGSTAB") 
	    default_fermi_inverter=BiConjugateGradientStabilizedInverter<fermi_field,gauge_field>;
	  else
	    error("unknown inverter");

	  sprintf(tmp,"%s_%f_%i%i.mdp", gauge_filename.c_str(),coeff_heavy["kappa_s"],b,j);
	  
	  if(is_file(tmp)) phi.load(tmp);
	  else {
	    psi=0;
	    if(on_which_process(lattice,0,0,0,0)==ME) { x.set(0); psi(x,b,j)=1; }
	    psi.update();
	    inversion_stats stats=
	      mul_invQ(phi,psi,U,coeff_light,
		       aprecision_heavy,
		       rprecision_heavy);
	    phi.save(tmp);
	  }
	  for(t=0; t<box[0]; t++) sum[t]=0;
	  forallsites(x) {
	    t=x(0); 
	    for(c=0; c<psi.nspin; c++)
	      for(k=0; k<psi.nc; k++)
		sum[t]=sum[t]+real(chi(x,c,k)*conj(phi(x,c,k)));
	  }
	  mpi.add(sum.address(),sum.size());
	  if(ME==0) {
	    fp=fopen("b_meson.log", "a");
	    for(t=0; t<box[0]; t++)
	      fprintf(fp,"%s %s %f %f %i %i %i %i %i %e\n",
		      gauge_filename.c_str(), smearing_type.c_str(),
		      coeff_light["kappa_s"],
		      coeff_heavy["kappa_s"],
		      a, b, i, j, t, sum[t]);
	    fclose(fp);
	  }
	}
  mpi.close_wormholes();
}
